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Abstract 

We study critical properties of the continuous Abelian sandpile model with anisotropies in 
toppling rules that produce ordered patterns on it. Also we consider the continuous directed 
sandpile model perturbed by a weak quenched randomness and study critical behavior of the 
model using perturbative conformal field theory and show the model has a new random fixed 
point. 

PACS: 05.65+b, 89.75.Da 

Keywords: Self- Organized Criticality, Sandpile models, Conformal field theory. 

1 Introduction 

The idea of self-organized criticality, introduced by Bak, Tang and Wisenberg [1], provides a useful 
framework for the study of non-equihbrium systems which dynamicahy evolve into a critical state 
without tuning of a control parameter. At critical state, these systems show scaling behaviors and 
this scaling behavior of the system is characterized by critical exponents 

The BTW sandpile model, renamed Abelian sandpile model after Dhar's work [3j, is the simplest 
lattice model that displays self-organized critical behavior. The Abelian structure of the model 
allows the theoretical determination of many of its properties |4j(5j. This model is usually defined 
on a square lattice. At each site of the lattice an integer height variable between 1 to 4 is assigned 
which represents the number of sand grains of that site. The evolution of the model at each time 
step is simple: a grain of sand is added to a random site. If the height of that site becomes greater 
than the critical height he = 4, the site will be unstable; it topples and four grains leave the site 
and each of the four neighbors gets one of the grains. As a result, some of the neighbors may 
become unstable and toppling continue. The process continues untill no unstable site remains and 
the avalanche ends. To achieve this, one should let some grains of sand leave the system, and 
this happens at the boundary sites. Every avalanche can be represented as a sequence of waves of 
the topplings such that each site at a wave topples only once [6]. While the scaling behavior of 
avalanches is complex and usually not governed by simple scaling laws, it has been shown that the 
probability distributions for waves display clear power-law asymptotic behavior [7]. 
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The scaling exponents of the system shows httle dependence on parameters such as the number 
of neighbors, however if we make the topphng rule anisotropic, then new universality classes may 
emerge: selecting a particular transport direction in BTW model, Hwa and Kardar [8] defined an 
anisotropic sandpile model such that the grains are allowed to leave the system only at one edge 
of the system. They determined the critical exponents with a dynamical renormalization group 
method. Dhar and Ramaswamy [9] defined a directed version of the BTW model and determined 
the critical exponents and the two-point correlation functions exactly in any dimensions. In |10| . 
the effect of anisotropy in a continuous version of sandpile model (Zhang model) is investigated. 
In this paper, a d-dimensional lattice is considered. This d-dimensional space is divided to two a- 
dimensional and {d — a)-dimensional subspaces. It is assumed that the energy (sand) is propagated 
differently for the two subspaces, but inside the subspaces the propagation of energy is isotropic. It is 
then shown that the peaked energy distribution and critical exponents of the distribution avalanche 
sizes are affected by the anisotropy. In jll] two variations of continuous abelian sandpile model are 
introduced, the directed model and the elliptical model. It is shown that the elliptical anisotropy 
does not change the universality class of the isotropic model whereas the critical exponents are 
sensitive to the directed anisotropy. Karmakar showed that in a quenched disorder sandpile model, 
the symmetric or asymmetric flows of sands in each bond determines the universality class of the 
undirected model ^12]. Also, a quenched disorder directed sandpile model has the same critical 
exponents with the BTW model when the local flow balance exist between inflow and outflow of 
sands at a site. Otherwise the model falls in the universality class of the Manna sandpile model 

The original isotope model could be represented with a conformal fleld theory known as c = —2 
theory [14j . When we insert anisotropy in toppling rules the rotational symmetry of the lattice is 
broken and the fleld theory associated with the model could not be conformal fleld theory. However, 
it may be possible to restore the rotational symmetry in large scale or statistically in an anisotropic 
sandpile model. To do this, one can introduce models that the toppling rules have some patterns 
on the lattice in a way that in larger scales there will be no preferred directions; that is, locally 
you have preferred directions which differ site to site in a regular pattern such that on larger scales 
the system look isotrope. Another possibility is to assume a quenched randomness for anisotropy 
in toppling of lattice sites; that is, we add anisotropy to the toppling rule of each site, however the 
amount of anisotropy and the preferred direction of anisotropy differs site to site randomly. In this 
way there may be no preferred direction statistically. 

The question we address in this paper is whether the universality class of these modifled models 
is different from the original sandpile model or not. We show that in some patterned sandpile 
models the universality class is the same as the isotropic Abelian sandpile model's universality 
class. However it turns out that the presence of disorder in a sandpile model may change the 
universality class of the system. This is done exploiting the replica technique; we consider this 
anisotropy as a perturbation to the original conformal field theory and use renormalization group 
to describe the perturbative behaviors of the system [15l [16] 

The plan of the paper is as follows: in next section we insert some anisotropics in the redistribu- 
tion of sands such that create some ordered patterns. We obtain the free energy function for these 
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models with using one to one correspondence between the recurrent configurations of ASM and the 
spanning tree configurations on the same lattice [T7]. The effect of these types of anisotropies on 
the critical behaviors of the system Theoretically and numerically is investigated. Next we consider 
a position dependent randomness in the toppling rule. Our procedure is based on the perturbative 
renormalization group approach around the conformal field theory describing the isotropic model 
and obtain the renormalization group equations for coupling constants. 

2 Patterned Continuous Sandpile Models 

It is known that the universality class of directed sandpile model is different from the ordinary 
ASM's [Hi [To]. In the directed model, the sand grains are always drifted toward preferred direction, 
say up-right corner. We would like to see if the directedness is introduced to the model only in 
small scales, is the universality class changed or not. To this end we add the directedness locally 
in a way that on average there will be no preferred direction towards which the sand grains move. 

Consider the continuous ASM on a square lattice composed of N lattice sites [?, ?]. To each 
site, a continuous height variable in the [0,4) interval is assigned. We divide the sites into two 
groups A and B, such that neighbors of one site in group A belong to group B and vice versa. We 
impose an anisotropic toppling rules for points of these two sublattice differently: when a toppling 
occurs in an ^-site 1 + e amount of sand is transferred to each of the right and up neighbors and 
1 — e amount to the down and left neighbor sites. In the case that a B site topples, 1 + e amount 
of sand is given to each of the left and down sites and 1 — e amount of sand is transfered to the 
right and up neighbors. Here, e is a positive real parameter less than 1 that controls the amount 
anisotropy. For e = we will have the isotropic model and e = 1 characterizes the fully anisotropic 
model. This toppling rule means that the A sites try to direct the avalanche towards up-left corner 
and the B sites try to direct the avalanche to down-right corner, thus on average the sands do not 
move in any specific direction. In Fig. 1 such a lattice is sketched. If a toppling occurs, the amount 
of sand transfered via thick lines is 1 -|- e and the amount of sand transfered via thin lines is 1 — e. 
It is clear that for e = 1 the sands are only allowed to move along one of the thick zigzag paths and 
therefore the system becomes essentially a set of one-dimensional sandpile models. The elements 
of the toppling matrix can be written in the following form: 




(1) 
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Figure 1: patterned ASM 
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Figure 2: Unit cells of the patterned ASM 



A first step to deduce the critical behavior of the system could be finding the free energy function. 
The closed form of the free energy is obtained by enumerating the corresponding spanning trees 
on the lattice. The formulation for enumerating spanning trees for general lattices is given in jl8j . 
We take the unit cells of two lattice sites as shown in Fig 2. Following the standard procedure, we 
obtain the free energy: 

f = ^ [ ^ dO I ^ d(P\ndetF{e,^) (3) 

OTT"^ Jo Jo 

where, 

F{e, (t>)= 41 - (a(0, 0) + a(l, 0)e*^ + a(-l, 0)6"^^ + a(0, l)e*'^ 

+a(0, -l)e-'^ + a(l, l)e^(^+'^) + a(-l, -l)e-*(^+<^)) (4) 

and a{n, n) are the 2x2 cell adjacency matrices describing the connectivity between sites of the 
unit cells n, n. 

a(-l,0)=a(-l,-l)=a^(l,0) = a^(l,l)= f ° ^~'] (5) 
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With a straightforward calculation one finds 

f = de /^^(i(/)ln(12-4e2-4(l-e2)cose-4(l+e2)cos(/)-2(l-e2)cos(^+0)-2(l-e2)cos(^-</))) 

OIT^ Jo Jo 

(6) 

Now it is seen that the model is equivalent to a free fermion 8- vertex model with weights {w{l), . . . w{8)} 
|19j that are related to e with the following relations: 

12-4e2 = w{l)^ + w{2f + w{3f + w{4f 

2(1 -e^) = wi2)w{4) - w{l)w{3) 

(l-e2) = w{5)wi6) - wi3)w{4) 

2(1 + e2) = w{2)wi3) - wil)w{4) 

= w{5)wi6) - wi7)w{8) 

The critical properties of the free fermion model are well known [19j . It is found that for all values of 
e the free energy function is analytical and the model shows no phase transition. It means although 
by inserting this kind of anisotropy some symmetries of the lattice are broken, but the broken 
symmetry operator is irrelevant and takes the system to the original critical fixed point. This fact 
can be checked by numerical simulations. We have simulated the model on a square lattice with 
sizes L = 64, 128, 256 and 512. After the system arrives at recurrent configurations, we began to 
collect data. At each size 10^ avalanches have been considered to derive the wave statistics. Fig 
3 displays the wave toppling distributions for different system sizes and three different values of e. 
A power law fit to these curves determines the critical exponent ri^^ defined as P^{s) ~ s~'^^ \ 
In Fig. 4 the extrapolated value of r for L ^ oo is obtained: t{oo) = 1.00 it 0.01 for e = 0.1, 
t(oo) = 0.99 ± 0.01 for e = 0.4 and t(oo) = 1.01 it 0.01 for e = 0.8. As we see, the wave exponents 
are independent of e and are consistent with the exact value of ri"'^ = 1 for e = [20] . 

It is possible to reformulate the partition function or the number of the spanning trees on the 
lattice in terms of fermionic path integrals. We place a two-component Grassmannian variable 
ipn = (V'l) V'2) on each unit cell n of the lattice. In this representation, the action of the field theory 
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Figure 3: Wave size distribution for e = 0.1, 0.4, 0.8 and for lattice sizes L = 64, 128, 256, 512. 
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Figure 4: The exponent t(L) is a linear function of 1/logL. The intersection with vertical axis 
gives Ts{oo) 

is written in the following form: 

2 

<n,n> i,j=l 

where a{n,n) are the adjacency matrices defined in ([5]). In the continuum limit, this action is 
obtained to be: 

r 2 

5 = / dxdy ^ [ii-ir+^i;lix)My) + 2e°^(l " e)d,i^iix)My) + 2e°^Vi(x)5,V'/3(y)] (8) 

o,/3=l 

where e"^ is the Levichivita antisymmetric tensor. At the first sight it may look strange that we 
have an action that has only first derivative in it, in contrast with the c = — 2 action that has 
second derivative terms. Even if we take the e ^ limit, it seems that the problem still exists. 
But if we look more closely, we will see that at least in the above limit one can write 'tp2 in terms 
of tpi and its derivative and then the second-derivative terms emerge. 

It may be argued that the above defined patterned system actually has a preferred direction; 
the zigzag paths join the down-left corner to up-right corner and not down-right to up- left. This is 
true, in fact the system has an elliptical anisotropy in the large scales and we know the elliptical 
anisotropy does not change the universality class [llj. It is possible to introduce other patterns 
in a way that the system be symmetric in the large scales. Fig. 5 shows such a pattern. In this 
model the thick lines characterize bonds that carry 1 + e amount of sand and the thin lines carry 
1 — e amount of sand after a site topples. Following the standard procedure, the free energy of this 
system is obtained: 

2 r2TT r2-!T p 

/ = ^ dO (i(^ln 132 - 136e2 + 4e^ + 2(l-e2)2(cos26l + cos2(/.) 

Ibvr^ Jo JO L 

-64(1 - e2)(cos <i) + cos 6) - 4(1 - €^f{cos{e + ^) + cos{e - (/)))] (9) 
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Figure 5: Order and symmetric patterned ASM 

which is again a smooth function and similar to the previous model and the self organized-criticality 
has the same universality class as of the undirected sandpile model. 

Up to now, we have observed that the patterns that do not produce a preferred direction in 
large scales do not change the universality class. In the next section we will consider a quenched 
random anisotropy to see if the universality class is changed or not. 

3 Random directed continuous sandpile model 

In directed continuous sandpile model (DCSM) introduced in [llj, it has been assumed that after 
toppling of a site, 1 + e amount of sand move to left (and up) and 1 — e amount of sand move 
to right (and down); that is there exists a preferred direction for the transportation of sands. In 
other words, the rotational symmetry is broken in this model. In the continuum limit. It turns out 
that the action of the theory assigned to the directed model is the action of c = —2 conformal field 
theory perturbed by the relevant scaling fields cp = —2686 and cj) = —26d6. As these operators are 
relevant, they grow under renormalization and take the system to a new fixed point [11] . 

In DCSM, e determines the strength of anisotropy and is in the interval (—1,1). Positive e 
means that the sand grains are pushed to the up-left corner and negative e means that they are 
pushed to the down-right corner. In this model the value of e is considered to be uniform through 
out the lattice. However, we may assume a statistical distribution for e, such that it can take both 
positive and negative values on different sites. The assumption that the mean value of e vanish, 
means that there will be no preferred direction statistically and the rotational symmetry will be 
restored to the model. The question is if such a modification takes the system to a new universality 
class or not. The assumption of a weak randomness allows us to determine the critical behavior of 
the model based on the pertubative renormalization group technique. 

In the continuous limit, the action of perturbed theory is given as: 

S = So+ f e{z,z){^{z,z) + ^{z,z)) (10) 
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Where 5*0 is the action of c = — 2 logarithmic conformal field theory. One can obtain the effective 
action using the replica method; that is, we have to take average of e on copies of the system 
and then find its limit when ^ 0. We assume that the e{z) at different sites are independent 
and have a Gaussian distribution on each site with a standard deviation equal to qq: 

{e{zi)e{z2)) = 9oS{zi - Z2) (11) 
The effective action then is expressed as: 

N N 

S = Y1 -^0,0 + 90 JZ^'^^^^' Z)4>biz, Z) + (t)a{z, Z)(j)b{z, z) + (l)a{z, z)(pbiz, z)) (12) 
0=1 •'^ aj^b 

Although the coupling constants of the field operators (f)<^, 4>4' ^^'^ f-i'e the same, as we will see, 
they have different RG equations. Therefore we distinguish the coupling constants of these field 
operators and rewrite them as go4>(j), Qq^^ and Qq^^ respectively: 



X] {ddH^aiz, Z)4>b{z, Z) + 9o^^(t>a{z, z)(t>b{z, z) + 9oct,^(l>a{z , z)(t)b{z, zfj 

N 

^ Y.^ab{z) (13) 

where the second line is an abbreviation of the first line. It is easy to see that the coupling constants 
9 are dimensionless; that is, they are marginal. Therefore to see if they are marginally relevant or 
not, we have to expand the partition function to the second order of 9. If it is marginally relevant 
we would like to see if it grows to infinity or will introduce a new fixed point, this means that we 
have to consider at least up to third order of coupling constants: 



/ E ^-b{z) + ^ / ^«b(^l) E ^-d{z2) + ^ / E E ^cd(^2) E ^e/(^3) 

+ --- = 9H / E '^»(^' ^)(l^b{z, Z) +94,4, jY^ 4>a{z, Z)4>b{z, Z) +94,4 Y ^»(^' ^)Mz^ ^) (14) 



a^b "^a^b " ^ aj^b 

To proceed, we have to know the contraction of fields in different possible ways. The calculation is 
done using operator product expansion (OPE) relations of the perturbing operators: 

(j){zi,Zi)4i{z2,Z2) = 7 + d(p{z2,Z2) + 2T{Z2,Z2) + . . . (15) 

[Zl - Z2y 

Hzi,zMz2,Z2) = , ^_ +d^{z2,Z2) + 2f{z2,Z2) + ... (16) 

[Zl - Z2) 

A.I -\ 1 , <^(2;i,^l) 4>{Zl,Zi) 

(p[Zl,Zl)4>{Z2,Z2) = 1 H z — + ■■■ (17) 

\Zl—Z2r Zl — Z2 Z1—Z2 

Where T and T are the components of energy-momentum tensor. 

At each order we contract all the fields using the above OPE relations and only keep a pair 
of (f) or ^ fields. While doing the integrations we have to perform regularization. We do the 
regularization in cut-off scheme: we assume the distance between any pair of integration variables 
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9<t><t> 

Figure 6: The RG flow of tlie model with quenched randomness. 

is restricted to be between a, the lattice constant, and L, size of the lattice. Up to the third order, 
the renormalized couplings are obtained to be: 

9<t><t> = 90H + 2a(A^ - 2)5o0<^5o</,0 + 2a'^{N - 2)[go^^gl^^{5N - 9) + go^^go^p^^N - 7)] (18) 

= 9000 + MN - 2)g^-^-^gQ^^ + 2a\N - 2)[go^^gl^^{5N - 9) + 5o009o00(3A^ " 7)] (19) 

500 = 5000 + MN - 3)(5o005o00 + 5o^^) + 8a[5o'<^^ {{N -2){N-1) + 2{N - 3)^) + 

500050005000 (3(iV - 2) (AT - 1) + 2(iV - Sf)] (20) 

where a = 47rln-^ and by the symmetry reasons, g^fp = g^^. In the limit = 0, we obtain the 
/3-functions up to third order: 

/^ff00 = = 167r5<^0500 - 16vra(95005j^ + 7^0) (21) 

= «^ = 247r(ff^^ + 55^)-32vra(54 + 64,^5^^) (22) 

It is clear from above equations that these fields are marginally relevant, however the coefficients 
of the terms proportional g^ are negative; hence the renormalization flow takes the system to a 
fixed point at g^fj, = = 0, g^;^ = (See Fig. [6|). In the new random fixed point, the rotational 
symmetry of the lattice restored so it is expected that the system show critical behaviors different 
from the deterministic directed model. 

We can compare our results with what Pan et al.|13j have found. In the patterned case, the 
outflow and inflow of the sand were balanced and we found that the universality class is not 
changed in such cases. On the other hand in the model with quenched randomness, there is not 
such a balance hence it is expected that the random fixed point belong to another universality class 
such as the universality class of the directed Manna sandpile model. We say it may correspond to 
Manna model because in this model there is randomness in the toppling rule, and we say it may 
be, because in Manna model the randomness is annealed but in our model it is quenched. 
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4 Conclusions 



In this paper we studied the critical behavior of the continuous sandpile model with the some 
patterned anisotropics in topphng matrices. Using the correspondence with the spanning trees, 
we obtained the free energy function for theses models. Both theoretic analysis and numerical 
simulations for the probability distribution of waves indicate that the anisotropic models are in the 
same universality class of the continuous sandpile model. 

Also we investigated analytically the cfFeet of quenched randomness on the critical behavior of 
continuous directed sandpile model. Our calculations is based on the perturbed renormalization 
conformal field theory and replica technique. Up to the third order perturbation, we obtained the 
renormalization group equations for the coupling constants of the perturbing fields. We showed 
that the perturbing fields are relevant and take the system to the new fixed point. 
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